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We present a unified kinetic theory that describes the finite-temperature, non-equilibrium dy- 
namics of a Bose-Einstein condensed gas interacting with a thermal cloud. This theory includes 
binary interactions to second order in the interaction potential and reduces to a diagonal quantum 
Boltzmann equation for Bogoliubov quasiparticles. The Hartree-Fock-Bogoliubov interactions in- 
clude the pairing field and are expressed as many-body T matrices to second order. The interactions 
thus include the correct renormalized scattering physics. This renormalized theory is automatically 
gapless. Thus, the excited Bogoliubov modes are naturally orthogonal to the condensate ground 
state. 
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I. INTRODUCTION 

Finite-temperature theories of Bosc-Einstcin conden- 
sation (BEC) have been a very active field of study. 
The goal is a unified description of a dilute, atomic gas 
of bosons in a harmonic trap in terms of a condensate 
mean field interacting with a thermal cloud. The success 
of the zero-temperature Gross-Pitaevskii (GP) theory in 
describing BEC experiments spurred interest in effects 
not contained in this framework. Examples are damping 
of collective excitations and condensate growth through 
collisional redistributions of thermal atoms. This redis- 
tribution cannot be treated in theories that are still of 
first order in the binary interaction, such as, for exam- 
ple, Hartree-Fock-Bogoliubov (HFB) theories. 

We present a kinetic theory of second order in the inter- 
action, which is formulated in terms of Bogoliubov quasi- 
particles and contains collisional terms beyond the HFB 
approximation. The HFB interactions are expressed as 
many-body T matrices to second order in the binary po- 
tential and thus include the correct renormalized scat- 
tering physics. This theory thus contains no ultraviolet 
divergences and has a gapless energy spectrum. We ex- 
tend the papers of Walser et al. [1, 2], which are the basis 
of the present results, by these two essential aspects. 

We go beyond theoretical approaches that drop the 
anomalous pair matrix rh — (aa) in the Popov ap- 
proximation [3-6]. Monte Carlo simulations based on 
the semi-classical Zaremba-Nikuni-Griffin theory [3] show 
very good agreement for experimentally observed damp- 
ing rates and response frequencies [7, 8]. However, recent 
experiments [9] and their theoretical explanations [10, 11] 
have shown that the pairing field plays an important role 
in Bose gases with resonance interactions. 

Keeping the pair matrix rh in this theory would cause 
ultraviolet divergences if we replaced the bare interac- 
tion potentials with the s-wave scattering length a s in 
the contact-potential approximation. The vacuum part 
of the pairing field's self interaction, for example, will di- 
verge, because the delta-function potential contains un- 



physically high energy contributions. These divergences 
can be resolved by writing the interactions in terms 
of scattering T matrices, which subsume the divergent 
sums over intermediate scattering states and correctly 
reduce to the scattering length a s in the zero-energy and 
-momentum limit [12-16]. Using these scattering T ma- 
trices, we can consistently eliminate all divergences to 
second order in the interaction. 

We thus obtain a renormalized HFB operator, which 
has a gapless spectrum. The zero-energy eigenspace is 
spanned by the condensate, and excitations with non- 
vanishing energy are thus automatically orthogonal to 
the condensate. This is similar to the renormalized gap- 
less HFB equations proposed in Ref. [14]. Other ap- 
proaches have to explicitly project the excitations orthog- 
onal to the condensate [12, 17]. Using the condensate as 
the ground state of this adiabatic basis also simplifies the 
complicated representation of a condensate band in the 
Gardiner-Zoller master equation formulation [18]. 

The Monte Carlo simulations of Jackson and Zaremba 
[7, 8] show that considering dynamic population-ex- 
change between the condensate and the thermal cloud 
and within the thermal cloud leads to good agreement 
with the observed response spectra. We include these ki- 
netic effects, thus going beyond collisionless descriptions 
[19, 20]. 

Our presentation builds on the papers by Walser et al. 
[1, 2] and we will begin by summarizing the kinetic equa- 
tions derived in these papers in Sec. II. In Sec. Ill, we 
examine the first-order HFB propagator and rewrite it in 
terms of second-order T matrices, by including second- 
order energy shifts and adiabatically eliminating the pair- 
ing field rh. This shows that the theory is explicitly gap- 
less and renormalized. We then make use of the gapless- 
ness and write the kinetic equations in terms of Bogoli- 
ubov quasiparticles, which are orthogonal to the conden- 
sate by construction. The results in Sec. IV show that 
the complicated collision terms presented in the Walser 
et al. papers can be simplified dramatically by a basis 
transformation. Practical calculations of the quantum 
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Boltzmann equation then will require only diagonal ele- 
ments of the quasiparticle population matrix. 



II. SINGLE-PARTICLE KINETIC EQUATIONS 

We present the Walser et al. formulation of kinetic 
theory [2], which was originally derived using a statis- 
tical operator approach [1]. The authors of this pa- 
per also recently showed a derivation of the same equa- 
tions [21] from the Kadanoff-Baym [22, 23] theory of 
non-equilibrium Green functions. They used the gapless 
second-order Beliaev approximation [4, 24] and showed 
equivalence to the work of Walser et al. under the Markov 
approximation. Recently, another independent deriva- 
tion [25, 26] connected this approach to kinetic theories 
by Morgan and Proukakis. 

Neglecting three-body and higher interactions, we can 
describe the weakly interacting, dilute gas by the follow- 
ing Hamiltonian 

H = H^ 1,2 'a\,a 2 , + 4> 1 ' 2 ' 3 ' 4 ' al^a^, (1) 

where # (0)1 2 = \2'} denotes the matrix ele- 

ments of the interaction-free single-particle Hamiltonian 



£(0) = JB_ + y ext(x) 
2to 



(2) 



with external harmonic potential V cy ±. The bosonic cre- 
ation operator a\, creates a particle in the state |1'), 
where 1' stands for a complete set of quantum num- 
bers, which label a constant, single-particle energy ba- 
sis, such as, for example, harmonic oscillator states or 
eigenstates of the GP equation. We use the summation 
convention for these abbreviated indices and indicate the 
single-particle basis with primes. 

The two-particle matrix elements of the binary inter- 
action potential Vbi n (x 1 ,x 2 ) are defined by 

1 ' 2 ' 3 ' 4 ' = \J dx i dx 2 (l'|x 1 )(2'|x 2 )t/ bin (x 1 ,x 2 ) 

x{(x 1 |3')(x 2 |4') + (x 1 |4')(x 2 |3')}.(3) 

These matrix elements are symmetric in the first and last 
two indices: 



il'2'3'4' _ j2'1'3'4' _ jl'2'4'3' 



(4) 



To determine the measurable quantities we want to 
calculate in this theory, we first define the condensate 
mean field a as the expectation value of the destruction 
operator 



a = a v |1') = (a v ) |1') . 
The total density matrix / is defined by 



(5) 



Subtracting the condensate density matrix 

f c = a®a* = a* 2 ,a v \l')®{2'\, (7) 

we obtain the density matrix of thermal atoms f = f—f c . 
The anomalous average to is split analogously 



m = {a v a 2 ,) |l') <g> \2') = m c + to, 



(8) 



in the condensate part m c = a <g> a and fluctuations to. 

This set of variables contains all possible combinations 
of up to two field operators. The course-grained sta- 
tistical operator parametrized by these variables is thus 
Gaussian, and we can use Wick's theorem to truncate the 
coupling to higher-order correlation functions, i.e., ex- 
pectation values of more then two field operators. This 
approximation is valid because of the diluteness of the 
condensed gas. In a dilute gas the duration of a collision 
event is short compared to the essentially free evolution 
between collisions, which allows higher-order correlations 
to damp. 

The procedure followed by Walser et al. [1] is then to 
write the Heisenberg equations of motion for the variables 
above and expand the expectation values using Wick's 
theorem. Walser et al. thus obtain the equations of mo- 
tion given in the following two sections [2] . 



A. Mean- field equations 

Since the anomalous fluctuations to couple the mean 
field a to its conjugate a* — a\, (l'|, it is convenient to 
write the generalized Gross-Pitaevskii (GP) equation in 
a two-by-two matrix form 



_ x = (-ffl + T<-T>) X , 
where the two-component state vector 



(9) 



* 1 „• I (10) 



is defined in terms of a = a v |1') and also contains the 
time-reversed mean field a* . 

The generalized GP propagator representing the re- 
versible evolution of \ is defined as 



n 



(ii) 



This symplectic propagator consists of the normal Her- 
mitian Hamiltonian 



= + u r + 2u f - & 



(12) 



where \i removes rapid oscillations of the mean field, and 
the symmetric anomalous coupling 



/=(a 2 ,a 1 ,)|l')®(2'|. 



(6) 



^ = Vfh. 



(13) 
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The energy shifts due to both the mean field and the 
normal fluctuations are given by the matrices 



[/ / = 2 1 ' 2 ' 3 ' 4 7 3 '2'|l / }®<4 / |, 



(14) 



whereas the first-order anomalous coupling-strength is 
given by 



Vrr, 



il'2'3'4' 



m 3 , 4 , |1') <g> |2') . 



(15) 



The second-order irreversible evolution, consisting of 
damping rates and energy shifts, is given by the collision 
operator 



T< 



(16) 



and its time-reversed counterpart T > 



-*iT<Vi, 



where o\ is the first Pauli matrix exchanging the positive- 
and negative-energy components of vectors. The matrix 
elements of the collision operator are given by 



T< 
T> 



= r /7(l+/) +2r /mm" ( 17a ) 

= r (l+/)(l+/)/ + 2r (l+/)mm" ( 17b ) 

= r mrrim* + 2r /m(l+/)' ( 17c ) 

= r mmm* + 2r (l+/)m/' ( 17d ) 



which is defined in terms of individual collisions T. These 
elementary collision processes arc defined explicitly as 



r /// 


_ o ,1'2'3'4' ,l"2"3"4" f f f M 
— of V 5 !) 73'l"/4'2"/4"2' I 1 


>®<3"|, 




C^l'2'3'4' ,1"2"3"4" / p 1 
= 80 0^ J 3 '2" m 4'3"/4"2' 1 


L'>®lO. 




o /1'2'3'4' /l"2"3"4" /• * 

= 80 (p n J 3 'i" m 4'4" m 2"2' 


l'>®<3"|, 


r 


1 'q'q'^_' 1 ft ntt ntt aH 

= 80 0,, m 3 , 4 „m 4 ,3„m 2 «2' 


|1')®|1"> 
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FIG. 1: The diagrams corresponding to the second-order 
terms T < in the GP equation (9). The dashed potential 
lines correspond to the symmetrized binary potential in the 
single-particle energy basis. The directed propagators repre- 
sent the normal density /, the remaining ones the anomalous 
average rh and its conjugate. Note that all diagrams are topo- 
logically equivalent, and only propagators are exchanged. 



sizable contribution. Note that the papers [1, 2] contain 
a sign error in the definition of Ae. For small i] we obtain 

1 -— — ► TrSJAe) + iV ri ^, (20) 

where V indicates that the Cauchy principal value has 
to be taken upon integration. The parameter r\ thus rep- 
resents off-the-energy-shell propagation after a collision. 
Most off-the-energy-shell coherences decay during subse- 
quent propagation, but, due to the finite time between 
collisions r c , energy cannot be conserved exactly, because 
7] has to be larger than the collision rate l/r c . 



Equations for normal densities and anomalous 
fluctuations 



The in-rates of the collision operator T are depicted in 
Fig. 1. 

In this theory, collisional interactions are considered 
to second order. The effect of higher order terms, which 
lead to a finite duration of a collision, can be modeled by 
introducing a parameter rj, such that every second-order 
collision operator contains dispersive as well as dissipa- 
tive parts from the complex-valued matrix element 



,1"2"3"4" 



,1"2"3"4" 



1 



rj — iAe 



(19) 



where the energy difference Ae — — y^ ir , -r c 2 „ } -r c 3 „ -r c 4 
has to be smaller than the energy uncertainty r\ to get a 



{el„+e%,) + e%,+el, 



The equations of motion for the fluctuation densities 
/ and rh are coupled and can also conveniently be writ- 
ten in terms of two-by-two matrices. To achieve this, 
we define the generalized single-time fluctuation-density 
matrix G < as 



G < 



f rh 
rh* (1 + /)* 



(21) 



where / = jW |1') ® (2'| and rh = m V2 > |1') ® |2') arc 
the matrix representations of the master variables. In 
Ref. [21], we use the property that this density matrix is 
the single-time limit of the corresponding time-ordered 
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two-time Green's function. This showed that the other 
time ordering is given by 



G ; 



(1 + /) m 
m* /* 



= a 1 G < *a 1 =a 3 + G < . (22) 



Here, we use the third Pauli matrix a 3 = diag(l, — 1). 
Note that our naming of the fluctuation-density matrices 
G < and G > is consistent with the two-time formalism in 
Ref. [21], but disagrees with Ref. [2]. 

The generalized Boltzmann equation of motion for this 
fluctuation-density matrix can be written as 

j t G< =-iSG < +r < G > -r > G < +H.c. (23) 

This equation has to be solved under the constraints 

a*f = and a*m = 0, (24) 

which force the fluctuations to be orthogonal to the con- 
densate. 

Again, the equation of motion (23) has two parts: 
The reversible evolution is governed by the Hartree-Fock- 
Bogoliubov self-energy operator 



EjV ^A 



(25) 



which in turn consists of the Hcrmitian Hamiltonian 

Sa/- = H {0) + 2U f c + 2U f - - (26) 



and the symmetric anomalous coupling 
£U = V m c + Vfh- 



(27) 



The irreversible evolution introduced by second-order 
collisional contributions now consists of the collisional 
operator 



r< = 



(28) 



and its time-reversed counter part T > = — uiT^ni. The 
diagonal components of the collision operator are defined 
as 



r A/- - r (/+/=)/(i+/) + r // c (i+/) + r /77 c ( 29 ) 

(f+f c )mm* + ^ fm c m* + ^ fmm c * } ' 
T M = r (l+/+/-)(l+/)/ + r (l+/)/<=/ +r (l+/)(l+/)/<= ( 30 ) 



+ \l+f+f)mfn* +^l+/)in'ra* + ^(l+/)mm< 

and the off-diagonal, anomalous components as 



r< — r + r _ + r ~ 

A (fh+m c )mrh* ' mm c rh* ' rhrhm c 



1 A 



(31) 

+ 2 { T U+f)rh{l+f) + r /m<=(l + /) + r /m/=}' 

r + r + r (32) 

(rh+m c )nifh* ' rhm c fh* ' fhfhm c * V J 

+ 2 { r (l + /+/<=)m/ + r (l+/)m<=/ + r (l+/)m/<= } ' 



FIG. 2: The diagrams corresponding to 
the second-order terms T < in the general- 
ized Boltzmann equation (23) . The dashed 
lines depict the symmetrized binary po- 
tential <f> in the single-particle energy ba- 
sis. The directed propagators represent the 
normal density /, the remaining ones the 
anomalous average m and its conjugate. 
The first column of diagrams is identical 
to those depicted in Fig. 1. The remaining 
diagrams each have one of the three prop- 
agators replaced with an open condensate 
line. 
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Both the diagonal and off-diagonal incoming rates are 
depicted diagrammatically in Fig. 2. For every term 
that appears in the collisional terms of the generalized 



GP equation in T< [Fig. (1)], we here [Fig. (2)] have 
three additional terms, where in each of them, one of the 
three fluctuating contributions is replaced with the cor- 
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responding mean-field quantity. This replacement rule 
can be seen in the Beliaev collisional self energies pre- 
sented in Ref. [21] and is a consequence of the fact that 
the Boltzmann equation (23) can be generated from the 
GP equation (9) by functional differentiation [27]. 

When the collision operator r< is multiplied by CP as 
in Eq. (23), we get terms like 



r /7(i+/)( 1 + ft 



(!+/)(!+/)/ 



(33) 



where the second part comes from the time-reversed con- 
tribution T > G < . The diagonal parts are exactly the in 
and out terms of the quantum Boltzmann equation for 
the single-particle distribution function /. The remain- 
ing second-order contributions couple to the anomalous 
fluctuations rh and do not have an analogue in the quan- 
tum Boltzmann equation. In Sec. IV, we will rewrite the 
kinetic equations presented here in terms of Bogoliubov 
quasiparticles. Then all collisional contribution take the 
form of Boltzmann terms. 



III. GAPLESSNESS 



T MATRICES 



Our goal in this section is to explicitly show that 
the kinetic equations (23) and (9) are gapless. This 
should on one hand be obvious, because our previous 
paper [21] showed them to be equivalent to the Kadanoff- 
Baym equations [22, 23] in the gapless Beliaev approxi- 
mation [4, 24]. The first-order HFB propagator S, which 
appears in the kinetic equations, is, on the other hand, 
known to exhibit a non-physical energy gap in the long- 
wavelength, homogeneous limit [28]. 

We resolve this discrepancy by including second-order 
collisional energy shifts "P{r} into the HFB propagator 
and adiabatically eliminating the anomalous average rh 
in the first-order anomalous potential V m (15). This up- 
grades the bare interaction potentials in the first-order 
propagators X and II to the real parts of many-body 
T matrices. We can then approximate these T matri- 
ces with the s-wave scattering length without incurring 
ultra-violet divergences. 

The upgraded HFB propagator S, where all binary 
interactions are written as many-body T matrices, is ex- 
plicitly gapless and thus obeys the Hugenholtz-Pines the- 
orem [29]. The propagator thus has zero-energy modes, 
which are completely specified by the value of the conden- 
sate a. If we use the non-zero energy Bogoliubov modes 
of £ as a basis for the thermal excitations, the excitations 
will automatically be orthogonal to the condensate. We 
follow this idea in Sec. IV. 



pairing field rh. We integrate the first-order equation of 
motion for the anomalous average rh 



d _ 



dt 



rh = -iY,j^rh - irh'Ejv - iY, A (l + /)* - ifT, A , (34) 



which is obtained by taking the fh component of the 
generalized Boltzmann equation (23) and dropping the 
second-order terms, because we want to substitute the 
result for rh into the anomalous potential V m and only 
keep terms up to second order. In stationarity, we solve 
for rh in the eigenbasis of Sa/S 



IW|eS > > = (e?,-/*)|eS,> 1 



and obtain 



1'2' 



Vl2 



(1 + /) 



, f y2"2' 
2'2" J~ J 1'2 //Zj .4 



2 M -(e° +e°,) 



(35) 



(36) 



as an adiabatic solution. Adiabatic here means that this 
solution only includes time- variations with characteristic 
times long compared to the duration of a collision. We 
use the Cauchy principal value V to indicate omission of 
the divergent term in an energy integral or sum. This di- 
vergent 8- function term gives rise to the imaginary part. 
We plug this result into the off-diagonal potential (15), 



V- 



1'2' 



tl'2'3"4'- 



4V- 



(1 + 2/) 4 „ 2 „ 



A3"2"3'4' 



2li - (e§„ + e°„) 



-m 



3'4'i 



(37) 



where we dropped the recursive V m term in the anoma- 
lous coupling E.4 in order to keep Eq. (37) at second 
order. We will discuss the recursive term in Sec. HID. 

We then recognize that we can write the off-diagonal 
element of the HFB propagator XU as the real part of a 
many-body T matrix 

£4 - V m c+V m = T m c(2fi), (38) 

which is defined by 

, ,,/ , rti 1 ' 2 ' 3 " 4 " (1 + 2 f) ^3"2"3'4' 

T 1 2 3 4 (e) = 20 1 2 3 4 +47>— -Jh"2" ( > 



e-(e°„+e°„) 

(39) 

The energies e° are dressed by the normal and mean- 
field shifts, but are not the full quasiparticle energies, 
because they do not include the effect of the pairing field. 
Contractions of this T matrix with anomalous averages 
are defined by 



T^'(e)=T 1 ' 2 ' 3 ' 4 '(e)m 3 , 4 ,. 



(40) 



The T matrix defined in Eq. (39) is a function of en- 
ergy through its last two indices in the sense that its 
argument e = e% + e\, . 



A. Off-diagonal potentials 



B. Diagonal Potentials 



Here, we update the off-diagonal potentials Vr m c +m \ in 
the HFB propagator S by adiabatically eliminating the 



In this Section, we want to redefine the diagonal po- 
tentials Ufc and Uf as the real parts of T matrices by 
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using the second-order energy shifts "P{r}. With P{r} 
we here denote the principal- value part of the collisional 
terms in Eqs. (28) and (16) according to Eq. (20). We 
will begin by considering the condensate potential. 

-i2U J * + p{2r /ef(1+f) + r fffe 

_2r / c (i+/)/ ~ r (i+/)(i+/)/ c } 

= -i2U f c-V{T 1{1+2f - )fc } (41) 

We here assume real eigenfunctions for the single-particle 
energy basis and do not include any F terms involving 
the anomalous average rh, because they are at least of 
order V^ 3 in /(Ae) 2 according to Eq. (36). The second-order 
terms that contain only normal fluctuations / are used 
in Eq. (46) to rewrite the fluctuation potential Uj. The 
term in Eq. (41) can again be written in terms of a many- 
body T matrix 

^ + |nr i(1+2/)/ c}=T / c, (42) 

which is given by 

(43) 

The slight difference compared to Eq. (39) is resolved 
when we assume diagonal quasiparticle populations 
P\2 = Pi ^Ii as wm be justified in Sec. IV: 

f 1 , 2 ,=UI,P [ U%r=f 2 , 1 „ (44) 

where U is the transformation matrix to the quasiparticle 
basis. Alternatively we note that the T matrix is essen- 
tially constant for energy differences up to the duration 
of a collision. Contractions of the T matrix with normal 
averages are performed according to 

T l'4'= T l'2'3'4' (e +e )/3 , 2 ,. (45) 

We now consider the fluctuation potential Uf, again 
include the truly second-order energy shifts, and obtain 

-i 2Uj + P{^ff{i+f) - r (i+/)(i+/)/} 

= -i2U f --V{r i{1+f)f } = T' f , (46) 

where we get a different T matrix defined by 

,1'2'3"4"/-1 i f\ J,3"2"3'4' 

T' im (e) = 2^' 2 ' 3 ' 4 '+47^ iVo "'ot 

(47) 

which does not have a factor of 2 in the intermediate- 
population term (1 + /). This difference is due to the 
fact that the mean field a is not bosonically enhanced. If 
we assume diagonal population fy 2 ' = Svi' fvv, which is 
not a good approximation in this basis, we reproduce the 
results of Ref. [30] for the GP equation to second order 
in the interaction potential. In particular, the factor of 
2 in their many-body T matrix in the term correspond- 
ing to Eq. (46) gets canceled with a negative term from 
adiabatically eliminating their triple average. 



C. Renormalized Propagators 

Using the T matrices defined in the previous sections, 
we can now rewrite the generalized GP propagator II 
given in Eq. (11) and the generalized Boltzmann propa- 
gator S given in Eq. (25). 

The Hamiltonian of the GP equation is now 

IIV = H (0) + T fe + 2T' f - n, (48) 

and the anomalous coupling 11^ vanishes, because of the 
identity 

V A a* = {T m c(2 M ) - V m c}a* = {I>.(2/i) - U f c}a. 

(49) 

This means that the coupling between a and a* vanishes 
and they are no longer independent quantities, and that 
we can without loss of generality treat a as real. 

When we write the Boltzmann propagator in terms of 
T matrices, we can explicitly show that the energy spec- 
trum is gapless; this theory thus fulfills the Hugcnholtz- 
Pines theorem [29]. The diagonal part of the propaga- 
tor S is now 

SV = H (0) + 2T fc + 2T' f - (j,, (50) 

and the anomalous coupling is 

£^T m c(2 M ). (51) 

We now show that this renormalized Boltzmann prop- 
agator £' has a zero-energy eigenvector, i. e., its spectrum 
is gapless. We begin by writing the generalized GP equa- 
tion for the condensate ground state 

{tiW + 7>(2 M ) + 2T£(2 M )} a = pa. (52) 

This equation can be written in terms of the renormalized 
GP Hamiltonian (48) as U_\fa = 0, where the energies are 
now measured relative to the adiabatic chemical poten- 
tial /i. An immediate consequence of Eq. (52) is that the 
quasiparticle ground state 

P a - C , (53) 

with a normalization constant C, is a zero-energy eigen- 
vector of the renormalized £', because of the identity 

T m c(2 M )a* =T f c(2fx)a. (54) 

These zero-energy eigenvectors of the HFB propaga- 
tor are proportional to the condensate mean field a. All 
non-zero-energy eigenvectors We^o 01 S arc thus au- 
tomatically orthogonal to the condensate, and we can 
use the complete set {a, We^o} as a basis to describe 
the condensate interacting with thermal excitations [31]. 
Other approaches to finite-temperature theories [12, 17] 
have to explicitly orthogonalize their bases using projec- 
tion operators. We discuss this new basis in Sec. IV. 
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D. Ladder approximation 



T l'2'3'4' (c) 



24> 



1'2'3'4' 



+4P 



fluctuations are by definition orthogonal to the conden- 
sate. Another motivation for transforming to a diagonal 
first-order Hamiltonian can be found in the numerical 
results of Walser et at; they find that the reversible first- 
order evolution leaves the quasiparticle populations con- 
stant. Thus, in the quasiparticle basis, only the second- 
order collisional terms change the populations. A more 
detailed account of the transformation to the quasiparti- 
cle basis can be found in Ref . [35] . 
0i'2'3"4"^j _^ 2/) 4 „ 2 „T 3 " 2 " 3 ' 4 '(e) Since the quasiparticles consist of the eigenvectors of 

the self-energy propagator we consider the propaga- 
tor's 2n by In (n is the number of single-particle states 
considered) eigenvector matrix W defined by 



We can extend the second-order T matrices introduced 
in the previous sections to include ladder diagrams to 
all orders by using consistency arguments. We first note 
that by keeping the recursive Vm term on the right side of 
Eq. (37), we can extend our definition of the off-diagonal 
T matrix to 



(4, +ei„) 



(55) 



This is the real part of the many-body T matrix in the 
ladder approximation. 

In the previous section, we showed that the HFB prop- 
agator £ is gapless with the T matrices defined to sec- 
ond order. If we used the ladder T in Eq. (55) on the 
off-diagonal while keeping the second-order T in Eq. (43) 
in the diagonal, we would find an energy gap of third 
order, because the cancellation in Eq. (54) only works if 
the two Ts are identical. However, the Hugenholtz-Pines 
theorem tells us that the full theory should again be gap- 
less. We thus conclude that we have to upgrade the T 
matrix on the diagonal of £ to the ladder approximation 
as well. 

We would like to finish our discussion of the scatter- 
ing matrices in terms of the single-particle energy basis 
with two remarks. First, as the Liouville-space formula- 
tion [32] of density-matrix evolution shows, the scattering 
should really be formulated in terms of Liouville-space 
scattering T matrices [33], which can be expressed in 
terms of Hilbert-space T matrices as 



T = TOl + l®T t + T(g)T t . 



(56) 



Since we only consider Hilbert-space T matrices, we thus 
would miss higher-order terms of the type TcgT^ , even if 
we included the full many-body Hilbert-space scattering 
matrix. The imaginary part of T gives rise to the inelastic 
rates in the kinetic equations in Sec. IV. 

Second, since the asymptotic states in this scattering 
problem are typically trapped harmonic-oscillator states 
for dilute, trapped atomic gases, we are strictly speaking 
dealing with bound-state R matrices [34] instead of T 
matrices. 



IV. QUASIPARTICLE KINETIC EQUATIONS 

A. Quasiparticle Basis 

We now want to write the kinetic equations (23) in the 
Bogoliubov quasiparticle basis. In this basis, the com- 
plicated and non-linear evolution due to the HFB propa- 
gator £' is replaced with a simple commutator with the 
eigenenergies and a slow basis rotation. As this theory is 
gapless, the E ^ quasiparticle states together with the 
condensate a form an orthogonal basis, i. e., the thermal 



E'TT = WE 



(57) 



at each time with the diagonal quasiparticle eigenen- 
ergy matrix E, which is labeled with the quasiparticle 
indices 1 = E\. The eigenvalue equations (57) are the 
Bogoliubov-de-Gennes equations. We decompose their 
solution W into two n by 2n matrices U and V* , 



W 



U 
V* 



u + U- 



(58) 



which are in turn split into n by n matrices for positive 
(it+, v*A and negative quasiparticle energies (it-, v*_). 
These quasiparticle eigenenergies 1 are the column in- 
dices and the original single-particle energies 1' the row 
indices, such that U — \U\^\^ V j^. 

Since the propagator 03E' is positive semi-definite, the 
eigenvalues E in Eq. (57) are real and come in positive 
and negative pairs of equal magnitude [36] . We argued in 
the previous section that there exists a zero-energy eigen 
vector P a defined in Eq. (53). This means that the prop- 
agator £' is degenerate and has a two-dimensional null 
space. The eigenvalue equation (57), however, does not 
yield the second linearly independent zero-energy eigen- 
vector. Instead, the associated vector Q a is given by [36] 



Pa 

Yj ' Qa = * A^' 
with a positive constant M. In our case, we find 



Qo 



-iC 



(59) 



(60) 



up to a normalization constant C . In order to find a 
complete set of basis states, we now define the quadrature 
components of P a and Q a : 



W +0 = -L(P a+l Q a ) 



W 



V2C 
iQa) = V2C 



(61a) 
(61b) 



If we now substitute these two states W +0 and IV~° 
into the zero-energy columns of W, we can normalize all 
eigenvectors using the symplectic norm 



W^a 3 W = cr 3 - 



(62) 
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The negative-energy states thus have negative norm, and 
the positive-energy states positive norm. In particular, 
the zero-energy vector W +0 has positive norm and thus 
belongs to the positive-energy part of W. This normal- 
ization fixes the constant C for the zero-energy states 
to 



C = 



1 



V2ak 



(63) 



where we defined the number of condensate atoms N c . 
The completeness relation for the quasiparticle basis also 
has symplectic structure 



2n, 



(64) 



where I271 indicates the 2n-dimensional unit-matrix. The 
symplectic structure of the orthonormalization and com- 
pleteness relations Eqs. (62) and (64), i.e., the appear- 
ance of the matrix er 3 in these expressions, guarantees 
that the transformation to the quasiparticle basis is 
canonical [36] . Canonical here means that the new quasi- 
particle operators obey the boson commutation relations. 

To further examine the structure of the quasiparticle 
states W, we consider the following symmetry of the HFB 
propagator £' 



£' = — (7lE'*(7l, 



(65) 



which holds according to its definition in Eq. (25). This 
symmetry implies [2, IV. B.] the following relation for the 
quasiparticle states W: 



W 



u + U- 

vl v* 



u+ v+ 
vl ul 



(66) 



This relation in Eq. (66) shows that the positive- and 
negative-energy eigenvectors of £' are not independent 
but related by 



W 1 = a x W- 



1* 



(67) 



Plugging W from Eq. (66) into the full completeness re- 
lation Eq. (64), we find a completeness relation for the 
independent elements of W alone: 



N, 



■aa 



E 

1>0 



1, 



(68) 



This is the basis we want to use for the kinetic equa- 
tions. The explicit split in a condensate mode a and 
the orthogonal fluctuation modes w 1>0 and v 1>0 is sim- 
ilar to the formalism of Ref. [31]. While the condensate 
mode a evolves according to Eq. (9) with the updated 
propagator II', we have to find an evolution equation for 
the quasiparticle populations. 

We thus write the generalized single-time fluctuation- 
density matrix G < defined in Eq. (21) in terms of quasi- 
particles as 



G< = WPW^ 



UPW UPV t 
V*PW V*PV t 



(69) 



where P is the not necessarily diagonal 2n by In quasi- 
particle population matrix. The time-reversed density 
matrix G > transforms according to its definition Eq. (22). 
From the same equation, we can also deduce that the 
quasiparticle population matrix P is Hermitian and ful- 
fills a similar identity: 

0-3 + P = oiPVi =a 1 P T a 1 . (70) 

This identity implies that P can be written as 



P = 



p q 

q* (1+pY 



(71) 



in a structure similar to G < in Eq. (21). Since the classi- 
cal mean field a only undergoes stimulated emission, the 
ground-state factors P+o+o = P-0-0 do not contain en- 
hancement terms. Thus, the zero-energy components of 
03 in Eq. (70) and similar enhancement terms associated 
with a actually have to be zero. 

We also represent the condensate (10) by a population 
matrix P c in the quasiparticle basis 



XX 



t = WP C WK 



(72) 



In order to maintain orthogonality between the con- 
densate and the thermal excitations, we have to demand 
an analogue to the constraint (24) in the quasiparticle 
basis. We here have to distinguish two cases. 

In the first case, the system evolves slowly enough that 
the ground-state of the adiabatic basis given in Eq. (68) 
is the true condensate state a(t) = a (real). We can then 
explicitly give the condensate matrix as 



(73) 



with a vector = <5j ± , such that all non-zero-energy 
components of P c vanish. To implement the orthogonal- 
ity constraint, we then explicitly set the E — ±0 elements 
of the quasiparticle matrix P to zero: 



01 



10 



0, for all 1. 



(74) 



In the second case, the system evolves too fast for the 
quasiparticle basis to follow. Then, the condensate ma- 
trix contains components of non-zero energy. In this case, 
the more general orthogonality constraint 



(a(t)\P = (a(t)\l)P r2 (2|=0 



(75) 



has to be fulfilled. This constraint is particularly impor- 
tant when the condensate is coherently excited in linear 
response. For adiabatic evolution, Eq. (75) reduces to 
Eq. (74). In the next section, we will find the evolution 
equation for P. 



B. Kinetic Equations 

With these ingredients, we can now try to obtain an 
equation of motion for the occupation number matrix P 
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from the kinetic equation (23). We use Eq. (69) to sub- 
stitute the fluctuation-density matrices G and obtain 



d_ 



(76) 



P = -i[E,P]_ + {w-^P + R.c.} 

+ {w- 1 r < w(<T 3 + p)-w- 1 T > wp + ii.c.y 



We are now left with the task of transforming the col- 
lisional contributions to the quasiparticle basis. To this 
end, we define new two-particle matrix elements in the 
quasiparticle energy basis by 



$ 1234 = 1'2'3'4' V^V^ulul. 



(77) 



These quasiparticle matrix elements have the same sym- 
metries as the original ones given in Eq. (4). Further- 
more, they fulfill 



*1234 



-3-4-1-2 



(78) 



The careful examination of symmetries in Ref. [35] 
shows that the collision operator Eq. (28) reduces to 
terms containing the completely symmetric matrix ele- 
ments 



1 



^1234 =J_Y^ 
V3^ 



(1234) 



(79) 



which arc defined as a sum over all index permutations 7r. 

The resulting kinetic equations for the quasiparticle 
populations P are 



d_ 
~dt 



P = -iEP + B w P + Cpp + C aP + ll.c., (80) 



with the non-adiabatic basis-rotation term 

R W-l dW M/t dW 

B W = W L — = a 3W T a 3 — , 



(81) 



and the Boltzmann collision rates between fluctuations 
Cpp = C< P (a 3 + P)-C> p P (82) 

= r P p (CT3+P) ((T3 + p) - T( a3+ p ){a3+ p )P p 

and between fluctuations and the condensate 

C a p = 3Tpc P ( (J3+ p- ) (<j3 + P) - 'ATpc^^+p^pP. (83) 

These collision rates are defined in terms of the quasipar- 
ticle collision operator 

vsTppp = r- :M '.-r ' ; /' , , /' , ;! /'; ; |3> ® <3'| . 

(84) 

In this operator, P can stand for any one of the three 
possibilities P, (a 3 + P), or P c , as needed in Eqs. 
(82) and (83). The approximately energy conserving ma- 
trix element tp^ is explicitly given by 



where the quasiparticle energies can be positive or neg- 
ative depending on their index. In an on-shcll scatter- 
ing event, this delta function forces two of the indices 
to be positive and the remaining two to be negative. 
This has to be considered in interpreting the collision 
terms Cpp and C a p, because there all the sums run over 
positive and negative indices. The principal- value part 
in Eq. (20), which appears in the single-particle kinetic 
equations, is absorbed in the T matrices in Sec. Ill and 
thus does not appear anymore in the quasiparticle equa- 
tions. 



o 



O 



FIG. 3: The incoming collision rates for collisions between 
thermal atoms Cp P and between a thermal and a condensate 
atom C^ P . The tight potential lines are now totally symmet- 
ric according to Eq. (79). The propagator lines represent the 
quasiparticle propagator P, which contains both the anoma- 
lous average rh and the normal density /. Because of the total 
symmetry of the interaction line, the three distinct conden- 
sate diagrams on each row of Fig. 2 reduce to one diagram 
with a weight of 3. 



To complete the presentation of the quasiparticle ki- 
netic equations, we write the generalized Gross-Pitaevskii 
equation (9) with the updated GP propagator II' and the 
second-order collisional contributions expressed in terms 
of quasiparticle populations P and obtain 



dt 



X 



-m'x+W<T 3 (Tpp {a3+ P ) -r (a3+P){lT3+ P ) p)W 1 x 



In Sec. V, we write this equation for a(t) alone. This 
equation together with Eq. (80) gives a complete descrip- 
tion in terms of Bogoliubov quasiparticles of a condensate 
coupled to a thermal cloud at finite temperatures. 



C. Orders of magnitude 

We now want to discuss the orders of magnitude of 
several of the quantities of this theory. This will sug- 
gest some approximations to the full quasiparticle kinetic 
equations (80). 

Wc first consider the basis transformation W. The 
completeness relation for the quasiparticle basis Eq. (64) 
tells us that 



1'2'3'4' 



1> 



1'2'3'4' 



nd v {E v +E- 2 , +E- 3 , +E l ,), (85) 



W = 0(1). 



(87) 
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For example, for high-energy eigenfunctions, the effect 
of the condensate becomes small, and the quasiparticle 
transformation reduces to 



1 and v 



0, as E 



oo. 



(88) 



Now, considering the basis-rotation term defined in 
Eq. (81) we find 



D 



w 



W 



_ x dW _ 0(1) 



dt 



dt 



(89) 



because the time-scale for population changes, which 
change the quasiparticle basis W, is limited by the time 
between collisions dt > t c . In equilibrium, the pop- 
ulations are constant due to detailed balance of the 
in and out rates. Thus, the net collision rate T = 
Cpp + C' a p + H.c. gives a better estimate for popula- 
tion changes dt w This also confirms that Bw = 
in equilibrium, since dW/dt = 0. 

We now show that the stationary solutions P of the 
Boltzmann equation (80) are diagonal. Considering the 
stationary solution 4jP\2 = of Eq. (80) for an off- 
diagonal clement with 1 ^ 2, we obtain 



. B W P + Cpp + C aP + H.c. (T 
= 1 E—El = °{ A-e 



(90) 



This shows that the off-diagonal elements of the quasi- 
particle population are small compared to the diagonal 
ones, which are of the order of the number of particles N, 
and vanish at equilibrium. 



The collisional terms of the kinetic equations (80) cor- 
respond to the imaginary part of a Liouville-space T- 
matrix (56). Examining the argument of the (5^-function 
in Eq. (85), which defines 'J,,, in comparison with the 
ladder T in Eq. (55) shows that (when the energies are 
correlated with the appropriate elements of W, see [35]) 
two of the energies in the denominator must be positive 
and two negative, which, together with the fact that P 
is diagonal close to equilibrium, leads to terms in the ki- 
netic equation of the form PP(1 + P)(l + P), etc. The 
kinetic equation (80) then becomes the quantum Boltz- 
mann equation for quasiparticle populations. The equi- 
librium solution is therefore the expected Bosc-Einstein 
distribution for the quasiparticles, as the steady-state so- 
lution of Eqs. (80) and (86) shows [2, Sec. V]. The inter- 
action matrix elements <fi in Eq. (77) can also be upgraded 
to Ts. These equations in terms of T are now consistent 
with an impact approximation treatment (with elastic 
scattering not contributing when T <g> terms are con- 
sidered, cf Eq. (56)). 



E. Conservation laws 

Because of the non- vanishing pairing field, the trace 
of the quasiparticle density matrix P is in general not 
equal to the number of excited particles in the system. 
Hence, the number of quasi-particles is not conserved. 
Our single-particle kinetic equations (9) and (23) do, 
however, conserve the mean total number 



D. Quasiparticle T matrices 

The T matrix in Eq. (39) (or the ladder extension 
Eq. (55)) has been obtained appropriately to second or- 
der. The second-order term with the factor of 1 in 
Eq. (39) is the divergent part, which is renormalizcd in 
Eq. (55) when replaced by the T-matrix. The term con- 
taining 2/ of Eq. (39) is a convergent many-body second- 
order term. It is thus reasonable to replace the T of 
Eq. (55) with 



T 1 ' 2 ' 3 ' 4 '^) = T} B 2 ' [i ' 4 '(E+2n)+8V 



1'2'3'4' , 



1'2'3"4" f i,3"2"3'4 
U"2"<P 

E — (Ey + E4/1 ) 

(91) 

where we replaced the single-particle energies in the de- 
nominator by quasiparticle energies, because the differ- 
ence is of higher order. We here use the following two- 
body T matrix 



r l'2'3'4' 
f 2B 



(e) = 20 1 ' 2 ' 3 ' 4 ' + 4V 



/ ) l'2'3"4"03"4"3'4' 



(92) 



which is given in terms of single-particle energies e sp de- 
fined by 

^ (0) |e i ?) = (£ + y ° xt(x) ) |e i ?) = $ |e i ?> • (93) 



(N) = a^a + Tr{/} = iV c + Tr{/} = const; (94) 

this can be proven explicitly by plugging these kinetic 
equations into (N) and canceling terms. We thus adopt 
a self-consistent procedure for the quasiparticle kinetic 
equations Eqs. (80) and (86), by requiring 



(N) = N C + Trjt/PC/t} = const. 



(95) 



This equation self-consistently constrains the number 
of condensate atoms iV c : in equilibrium at tempera- 
ture the quasiparticle matrix P consists of q = 
and 



p(E) 



e 0{E-,j.) _ I 



(96) 



according to Eq. (71) with the chemical potential fi given 
by the GP equation (52). As temperature tends to zero, 
we obtain the usual corrections for the anomalous aver- 
age and condensate depletion [12]. Note that the number 
of excited atoms is not given by the trace of p, but has 
to include the basis transformation U as discussed above. 
If we drop the basis-rotation term Bw m numerical sim- 
ulations, we will incur number-non-conservation on the 
order of T, while away from equilibrium. 

Since we use a Markovian collision integral with a 
damping function of finite width 77 in order to include 
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off-the-energy-shell propagation, this theory is not ex- 
actly energy conserving. Markovian theories fail to track 
the decay of initial correlations and thus do not account 
for the decay of the correlation energy [37]. In our case, 
with a self-consistent 77, energy is conserved to order 77, 
which is consistent with the order of approximation. For 
a detailed discussion of these memory effects and how 
they affect the conservation laws see [38]. 

V. SUMMARY 

We summarize the main results of this paper given in 
Sees. Ill and IV. We formulate a kinetic theory in terms 
of Bogoliubov quasiparticle modes W, which are defined 
by the eigenvalue equation for the renormalized Hartree- 
Fock-Bogoliubov propagator S' 

( m+2T f ^ 2T r» T -;)w = WE. (57) 

The T matrices are defined to second order in Eqs. 
(39) and (47) and Eq. (55) gives an extension to the lad- 
der approximation. The Gross-Pitaevskii equation 

{AW + 7>(2 M ) + 27^(2 M )} a = m (52) 

for the ground state a is contained in Eq. (57), because 
the renormalized £' is gapless. The GP equation defines 
the value of the chemical potential \i. To find a com- 
plete basis W, we have to find the second zero-energy 
mode and form quadrature components as discussed in 
Sec. IV A. 

We find the following Boltzmann equation for the ther- 
mal excitations in terms of Bogoliubov quasiparticles: 

jP = -i[E,P] + {a 3 W^a 3 ^-P + K.c.} (80) 

+ { r pp( ( T3+p)( cr 3 + p)- r( CT3+ p)( CT3+ p)pP + H.c.j 

+ 3{r P cp (CT3+P) (a3 + P) - r pc(a3+P)P p + H.c.}. 

The basis- rotation term containing 4rW can be dropped 
for adiabatic evolution. The quasiparticle density ma- 
trix P is diagonal close to equilibrium, and its elements 
obey an equilibrium Bose-Einstein distribution as the 
second and third lines show. The collision terms contain- 
ing the general condensate matrix P c defined in Eq. (72) 
represent population exchange between the thermal cloud 
and the condensate. They are balanced in the following 
equation for the condensate 

ja{t) = {HW+T fc +2T' } -»}a{t) (86) 

+ Ua 3 {T P p^ 3+P) - r( CT 3 + p)( CT3+ p)p}a3 
x{tfa(t) -V T a*(t)}. 



In general, a(t) can be different from the a used as the 
ground state of the adiabatic basis. A non-adiabatic 
change in a driving force, for example, would cause a(t) 
to change quickly. 

The coupled Eqs. (80) and (86) have to be solved under 
the orthogonality constraints (75) or (74) depending on 
whether the evolution is adiabatic or not. Furthermore, 
the total particle number has to fulfill the self-consistent 
constraint Eq. (95). 



VI. CONCLUSIONS 

We have extended the non-equilibrium kinetic theory 
of Walser et al. [1, 2] in two important respects. First, 
we write the binary interactions as many-body T ma- 
trices to second order in the interaction by subsuming 
ultra-violet divergent terms. This procedure removes di- 
vergences caused by the anomalous average rh and con- 
tained in the second-order terms. We can then replace 
the low-energy limit of the T matrix with the s-wave 
scattering length. This is to our knowledge the only con- 
sistent treatment of these difficulties associated with the 
anomalous average in a theory that includes second-order 
collisional terms. The updated Hartree-Fock-Bogoliubov 
propagator £' is then gapless, which greatly facilitates a 
consistent treatment of the condensate coupled to ther- 
mal fluctuations. 

The second extension of the Walser et al. theory makes 
use of the gapless HFB propagator by using its quasi- 
particle eigenmodes to parameterize the thermal fluctu- 
ations, which are then automatically orthogonal to the 
condensate mean field. We find that this basis greatly 
simplifies the second-order collision terms of the Walser et 
al. theory. Another important result reported here is 
that, in equilibrium, the Bogoliubov modes diagonalize 
the quasiparticle population matrix P. This means that, 
close to equilibrium, the second-order terms can be eval- 
uated in n 4 operations, where n is the number of energy 
levels considered. This is a vast improvement compared 
to n 8 operations for a basis that docs not diagonalize 
the population matrix, such as, for example, the single- 
particle basis used by Walser et al. 
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